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Abstract 

This paper describes a technique to determine the linear well- 
posedness of a general class of vector elliptic problems that include 
a steady interface, to be determined as part of the problem, that sep- 
arates two subdomains. The interface satisfies mixed Dirichlet and 
Neumann conditions. We consider "2+2" models, meaning two in- 
dependent variables respectively on each sub domain. The governing 
equations are taken to be vector Laplacian, to be able to make ana- 
lytic progress. The interface conditions can be classified into four large 
categories, and we concentrate on the one with most physical interest. 
The well-posedness criteria in this case are particularly clear. In many 
physical cases, the movement of the interface in time-dependent situa- 
tions can be reduced to a normal motion proportional to the residual in 
one of the steady state interface conditions (the elliptic interior prob- 
lems and the other interface conditions are satisfied at each time). If 
only the steady state is of interest, one can consider using other resid- 
uals for the normal velocity. Our analysis can be extended to give 
insight into choosing residual velocities that have superior numerical 
properties. Hence, in the second part, we discuss an iterative method 
to solve free boundary problems. The advantages of the correctly cho- 
sen, non-physical residual velocities are demonstrated in a numerical 
example, based on a simplified model of two-phase flow with phase 
change in porous media. 
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1 Introduction 

Free boundary problems (FBPs) have motivated several studies in 
the past due to their vast applications in fluid flow, phase change 
models and other fields. Some classical free boundary problems are: 
the dam seepage problem [BJ, incompressible two-phase flow Q3] 
(i.e. falling droplets or rising bubbles [12], the Alt-Caffarelli problem 
[2], the classical Stefan problem [6j, etc. From a mathematical point of 
view, FBPs are boundary value problems with an unknown boundary. 
The motion (unsteady case) or position (steady case) of the boundary 
has to be determined together with the solution of the given partial 
differential equations on one (free surface) or both sides (interface) of 
the free boundary. The coupling of the free boundary to the interior 
is always nonlinear [13], and thus FBPs are often not easy to solve. 

The common structure of all the examples above is that at steady 
state, they all have second order elliptic governing equations for m 
unknowns on one side and n unknowns on the other side of the free 
boundary. We denote this situation as an "m+n" problem. There 
must then be m + n + 1 Dirichlet-Neumann conditions at the interface. 
For fourth order problems such as biharmonic equation, we would 
need 2m + 2n + 1 conditions at the interface to determine an "m+n" 
problem. By this generalization, many more complicated problems 
can be formulated. However, fourth order problems and second order 
elliptic systems other than vector Laplacian are not considered further 
in this present study. 

To solve free boundary problems numerically there are three main 
kinds of methods: capturing methods, front tracking methods and 
level set methods. Capturing methods are based on Eulerian for- 
mulation and the problems are reformulated and solved in the whole 
domain. The interface location is recovered from the discrete solution. 
In these methods, the interface conditions are not specified explicitly. 
A classical example is the Enthalpy method [6j. The alternative is to 
discretize the interface explicitly [16] [18] or via a level set approach 
|llj . In both these cases, the interface conditions (and interface ve- 
locity for time-dependent problems) are implemented explicitly. This 
can be done by considering the domains as disjoint and discretizing 
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the equations and interface conditions directly [8] or by discretizing 
the entire domain and modifying the discretization near the interface 
|10| . The latter approach combined with modern level set techniques 
can also be considered a capturing method. 

In a few cases, steady state solutions can be obtained directly using 
shape optimization [THJ, [§]. Another approach to reach steady state 
solutions is to solve the transient time-dependent problem to long 
time. In the class of problems we consider, the normal motion of the 
interface is driven by the difference of solutions on either side of the 
interface. For example, solidification boundaries are driven by a Stefan 
condition |6j. It should be made clear that we are not considering here 
the problem of dendritic growth [IJ |4j which in our framework would 
be ill-posed, but is regularized by higher order (curvature) terms. We 
consider only well-posed problems which do not need regularization. 
In the physical example of section [31 the interface is moved by the 
net mass flux at the interface. Those are typical conditions that give 
normal velocity of each point on the interface. In these examples and 
others of physical interest, the normal interface velocity is given by 
the residual in one of the steady state interface conditions. In this 
paper, we apply the residual velocity method (a tracking technique) 
f(ZJ: 

1. Given an initial interface T and tracking points {x : x € T}, 
solve numerically a fixed boundary value problem satisfying all 
of the steady interface conditions but one. 

2. Substitute the solutions into the unsatisfied condition and find 
the residual i?(x). 

3. If the residual is larger than a tolerance, explicitly evolve the 
interface using the residual as a normal velocity of the tracking 
points: 

V„ = i?(x) 

4. Repeat the process until the residual is less than a given toler- 
ance. 

At termination, all interface conditions are approximately satisfied. 
This is an example of a value method [15] (since it uses only val- 
ues of the solution not shape derivatives to update the interface). In 
this algorithm, the choice of residual velocity is not specified, so we 
can choose different interface conditions and their combinations as the 
residual velocity, not only the real, physical velocity. As a result, we 
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can show both analytically and numerically that the residual velocity 
chosen by our criteria has better numerical properties than the phys- 
ical one. Specifically, our method can use time steps independent of 
mesh size in the explicit step 3 above. Because of the use of an artifi- 
cial interface velocity, our method provides accurate solutions only at 
steady state. 

There are few publications addressing the analytic theory of the 
general class of the steady free boundary problems considered here. 
In [7J the authors present a linear theory of the well-posedness of the 
general class of "1+1" Laplacian problems. Our work is an extension 
of [7J. Specifically, the ideas considered there for the scalar case are 
extended to the "2+2" case. An analytic simplification to the well- 
posedness is identified in a physically important class of problems. 
The class is represented by a simplified model of two-phase flow with 
phase change in porous media. In addition, a different numerical im- 
plementation of the corresponding residual velocity method is given. 
While the numerical method is not as general as that of [7J it does 
show that accurate results can be obtained without grid refinement at 
the interface. Attaining this in more general codes is a goal for future 
work. 

The framework of this paper is following: In section we discuss 
the four classes of second order "2+2" problems with different com- 
binations of Dirichlet-Neumann conditions. A physically interesting 
example, a simplified model of two-phase flow with phase change in 
porous media, is presented in section El In section [U well-posedness 
criteria of "2+2" Laplacian equations from [7J are reviewed and ap- 
plied to the example and generalized. The extension to 3-D problem 
has been shown as well. Then in section [5j the example of vector 
Laplacian type free boundary problem are computed numerically, us- 
ing a finite difference method in mapped coordinates (a so-called ficti- 
tious domain method [9j). The performances of both physical velocity 
and our residual velocity are compared, and the results agree with the 
analytical predictions. 

2 Class Division of "2+2" Problems 

Consider a two-dimensional model, when m = n = 2. We have a 
"2+2" vector elliptic problems for variables u ± ,p ± , assuming that 
Li(u + ,p + ) = in the upper subdomain D + and \j-z{u~ ,p~) = in 
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Figure 1: The model of second order "2+2" problems. Li 2 are linear, second 
order elliptic operators and there are five interface conditions represented by 
G. 



the lower subdomain D~ , as shown in Fig. [TJ There are five conditions 
at the interface. Here, we consider Li and L 2 to be general, vector, 
linear, second order elliptic operator but in our analysis below we 
consider only the case of L12 both being the vector Laplacian. 

The five Dirichlet-Neumann boundary conditions at T in matrix 
form, will be: 

GU = b. (1) 

where G is a 5 x 8 matrix, and U T = (n„ , u~, p^, Pa: n+ "> u ~> P + i p ) 
is a vector composed of all unknown variables and their normal (n 
pointing from D~ to D + ) derivatives at the interface. The matrix G 
can be split naturally into two parts: 

G = (Gn I Gd) 

where Gn and Gd are 5x4 sub-matrices, representing respectively 
Neumann conditions and Dirichlet conditions. 

We use the rank r of Gn (or equivalently the number 5 — r of pure 
Dirichlet conditions) as the criteria to classify possible forms of G. As 
a result, there are four classes of boundary conditions corresponding to 
r=4, 3, 2, 1 below. Normal forms can be obtained using the following 
allowable operations: 

1. Row operations on G. 

2. Simultaneous swap of columns 1 and 2, 3 and 4, 5 and 6, 7 and 
8. This swap is accompanied by a sign change of columns 1-4. 
(corresponds to a relabelling of D~ and D + ). 

3. Column operations between columns 1 and 3 with identical op- 
erations on columns 5 and 7 (linear combination of unknowns in 
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£>+). 

4. Column operations between columns 2 and 4 with identical op- 
erations on columns 6 and 8 (linear combination of unknowns in 
D-). 

The following generic forms result: 

Class A: One pure Dirichlet condition, and thus the rank of Gn is 
4: 
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Class B: Two pure Dirichlet conditions, and the rank of Gn is 3. 
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where the last two rows are linearly independent but no further 
reduction is possible. 

Class C: Three pure Dirichlet conditions: 





/ Gi,i 


£•1,2 


£•1,3 


£■1,4 


| G h8 


\ 




G2,l 


G2,2 


£•2,3 


£•2,4 


| G 2 [ 8 




G = 














| 1 G 3 , 8 


















| 1 G 4 , 8 






V 











| 1 G 5 ' 8 


/ 



where the first two row vectors restricted to the left block are 
linearly independent but no further reduction is possible. 

Class D: Four pure Dirichlet conditions: 



G = 



/ 1 


£•1,3 


£•1,4 


\ 











10 











10 











10 


V 








1 / 



6 



Class C is arguably of most physical interest. If the vector el- 
liptic equations represent conservation of two quantities with fluxes 
proportional to gradients, then the conservation of these quantities at 
the interface are represented by exactly two pure (Gig and G^s both 
zero) Neumann conditions. We define the subclass C to be the one 
with these additional conditions. In the following section, we present 
a physical example in this class originally derived from [3]. Class D 
always leads to a single pure Neumann condition. We define sub- 
classes A and B to be the ones with four and three, respectively, pure 
Neumann conditions. 



3 Physical example 



We consider now a physical problem involving two-phase flow with 
phase change in porous media. Consider a closed system consisting of 
a sand pack, water vapour and liquid water. If the system is heated 
at the top, vapour is formed, and in some cases, two distinct zones 
with free boundary will be observed. In the upper zone D + , there is 
only vapour, while in the lower zone D~ , two-phase vapour and liquid 
water can be found. This model is studied in [3]. 

In the upper region D + , the variable P v (x,t) and temperature 
T + (x,t) describe the vapour region, and the governing equations are: 

AT+ = (2) 
V-(p v ii„) = (3) 

where p v and u v denote density and velocity of the vapour respectively, 
which are the functions of temperature T + and pressure P v . According 
to Darcy's law, we have: 

u,, = VP„ 

Pv 

with k the permeability and fi v the dynamic viscosity. If we assume 
the gas is ideal, then: 

__ M P v 

where M is molar mass of water, R is the universal gas constant. So 
the equation ([3|) becomes: 

V • (^rVP v ) = (4) 
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Then ([2]) and © give an elliptic problem for P v and T + . 

Flow in the lower two-phase region is described by temperature 
T~ and saturation s(x, y). The saturation is the fraction of pore space 
occupied by liquid water. Energy conservation and mass conservation 
give us: 

V • (KVT~) - h vap V ■ (p v u v ) = (5) 
V • (piui + p v u v ) = (6) 

where u;,u„ and p v are functions of T~ and s. Darcy's law includes 
the relative permeabilities as well as the capillary pressure, which we 
assume are functions of the s: 

u v = - — (l-s) 3 VP sat (T-) 

Pv 

Ul = --s 3 V(P sat (T-) - P c (s)) 
Ml 

where p v ,Pi are respectively dynamical viscosities of the vapour and 
liquid, P sa t is the saturation vapour pressure, and P c = P v — Pi is the 
capillary pressure. Since evaporation and condensation rates scale to 
be very large in a porous media, we assume P v = P sa t(T~). Therefore 
in the two-phase region, we have an elliptic problem for T~ and s. 

On the free boundary, There are five Dirichlet-Neumann conditions 
required for the problem of four variables: 



s = 
[T]=0 
[Pv] = 

(p v u v ) + • n = (p v u v + piUiY ■ n 
[KT n ] = h vap (piui) ■ n 



saturation is zero at the interface 
temperature is continuous 
vapour pressure is continuous 
mass conservation 
heat conservation 



where [•] denotes the difference between counterparts on each side of 
the interface. For example: [T] = T + — T~ . 

In order to make progress on this problem in our framework, con- 
siderable simplification is needed. Specifically, we need to simplify 
the problem to get a vector Laplacian operator on each subdomain 
and linear interface conditions. The model derived below lacks sev- 
eral aspects of physical interest and mathematical difficulty. However, 
it does retain the underlying structure of two elliptic systems coupled 
by five mixed Dirichlet-Neumann conditions at the interface. First, 
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we set P v /T + = 1 in , which means the vapour density doesn't 
change. This leads to Laplacian equation of pressure: AP = 0, then 
in the upper region D + , we have: 

AP = 0, AT+ = 

In the lower region D~, we simplify using the assumption that 
evaporation or condensation only happens at the interface. Moreover, 
we assume the relative permeability of vapour and liquid in D~ is 
independent of saturation s, the capillary pressure P c is linear in s 
and saturation vapour pressure P sa t is linear in T~. In the end, ([5]) 
and © can be reduced to two Laplacian equations: 

As = 0, AT' = 0. 

Correspondingly, the interface conditions are given: 



s 


= 


saturation is zero across interface 


[T] 


= 


temperature is continuous 


P 


= T~ 


Vapour pressure is continuous 


KT n ] 


= S n 


heat flux is continuous 


Pn 


= T n + s n 


mass conservation 



Note that in the example, the interface is featured by s = 0, and 
its movement is driven by the mass flux going through the interface, 
so we say the physical residual velocity comes from the last Neumann 
condition. The five conditions can be written in matrix form ([TJ where 
G and U are respectively: 
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and 

U T = (T+, T~ , P n , s n , T+, T- P, s) 

Clearly, the boundary conditions are composed of three pure Dirichlet 
and two pure Neumann conditions, which belongs to Class C identified 
in the last section. 
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4 Linear Well-Posedness and Residual 
Velocities 



In the first subsection below we review briefly the theory developed in 
[7] targeted to the "2+2" case. Later, we identify situations where the 
analysis gives particularly clear insight. These include the physically 
important case C identified in the previous section. Application of 
the theory to the example above is done, and a residual velocity with 
superior numerical properties is identified. 

4.1 General Theory 
4.1.1 Base Solutions 

We find base solutions in 2D (coordinates x and y) with a flat interface 
y = 0. Take Uq = (u + , u~ ,p + ,p~) to be the base solution of the FBP 
satisfying: 

Au = in y > and y < (8) 
with interface conditions: 

G Uo = b on the interface y = (9) 

where Un = ( U ° n ) • 

We consider base solutions independent of x with b constant, so 
(|8j) can be reduced to uo^j, = and uo is linearly dependent on y. 
That yields: 

u = y r 1 + r° 

where 

is a constant vector solving (|9|). The base solution is composed of a 
particular solution and homogeneous solution 

U = r H + r P , 

with Grn = 0. Hence, yh belongs to the nullspace of G, and M{G) 
has rank 3. 
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In order to specify a unique solution, we need three more conditions 
denoted by a vector q = (q\, q2, q?,) T , corresponding to the ordering of 
a basis of vectors of M(G), so 

Uo= ( ) =AT(G)-q 

where we mean (with a slight abuse of notation) by M(G) the 8x3 
matrix with a basis for this space in columns. 

Consider the base solution found above to describe the local be- 
haviour of the interface near a point, where the coordinate system has 
been rotated so that the local tangent plane is y = 0. This local base 
solution depends on the local data b and also the vector q, which we 
refer to as the global fluxes. This represents the dependence of the 
local solution on far field conditions. 

4.1.2 Linear well-posedness 

We consider how small perturbation of boundary conditions about a 
flat interface affect the solution: 

GU = b + efe iax (10) 

where f is a constant, a is the Fourier mode, and e is a small parameter. 
As a result, we will obtain a new solution and new free boundary. If 
we can show that the linearized solution depends continuously on the 
data f for each a we say the problem is linearly well-posed. 

After we perturb the interface condition using (|1U|) . we expect the 
new solution u and free boundary y = r/(x) to be in the form: 

u(x, y) = u (y) + eui(x, y) + 0(e 2 ) 
rj(x) = 6771 (a;) + 0(e 2 ) 

where uq, rfo(x) are the base solution from section [4.1.11 above. Fol- 
lowing [7J it is known that linear perturbations are of the form: 

m (x) = me iax (11) 
u(x,y) = ce Tlaly e iax (12) 

where (I12p arises from the separable solution of the Laplace opera- 
tor in half planes and the sign is chosen appropriately depending on 
whether the corresponding variable is in the upper or lower domain, 
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respectively. Again following [7J, if the forms (jlip and (|12p are put in 
to (jlOp the first order term results in the following algebraic condition: 



GMc + r}iGU 0n = f 



(13) 



where 
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u 0yy 
u 0y 



evaluated at y = 0. This can be written as: 

U n = cSUq 

where Uo is the vector from the base solution (pQ) and S is a shift 
matrix that moves rows 1-4 to rows 5-8 and then zeros the first four 
rows. 

Considering (fl~3|) . the well-posedness criteria can be reduced to the 
existence and uniqueness of (0,771) for every f and a. Consider the 
matrix GM. It is a 5 x 4 matrix with rank 4. We denote a vector in its 
left nullspace by w = JV[(GM) T ], Thus to ensure the well-posedness 
of the problem, we require: 



w T GU 0n + 



(14) 



This is the algebraic condition for linear well-posedness for the general 
class of steady free boundary value problems we have considered. 



4.1.3 Residual velocity choice 

In our numerical algorithm, we solve a fixed boundary problem each 
iteration with some of the boundary conditions in G satisfied, while 
the residual in the remaining condition is used as the normal velocity 
to evolve the interface. With the well-posedness criteria of the prob- 
lems established above, we can discuss the stability of the numerical 
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methods for the free boundary problem. We consider residual velocity 
methods as discussed in the introduction. The only time dependence 
is in the interface position y = r/(x,t) which in the linearized setting 
above reduces to 

y = fje iax e xt 

The stiffness A is a function of the choice of residual to use as velocity 
and also the wave number a. If for a given residual velocity choice, 
Re{\) < for all a, the resulting method will converge to the steady 
interface solution. We show how in addition residual velocities can 
be chosen in some cases so that A is bounded independently of a, 
which suggests that time steps for the resulting method can be chosen 
independently of the spatial discretization. This is verified in the 
example computation in the next section. 

Again, following [7] the linear analysis leads to 

A= wTG T °°" (15) 

where v denotes the linear combination of residuals used as the normal 
velocity (the orthogonal combinations are used as interface conditions 
at each time) and as before w is the left nullspace of GM. 

Note that in the numerator of (|15p is the well-posedness term (jl4j) . 
We assume we are only trying to compute solutions to well-posed prob- 
lems, so can assume this term does not change sign. Thus, the stability 
of the scheme depends only on the denominator which depends on a 
straight-forward way on the velocity choice v. 



4.1.4 Discussion 

The well-posedness condition (I14j) and stability condition (115j) are 
powerful analytic tools. The influence of the boundary conditions 
(through G), the elliptic problem (through w coming from M), and 
boundary and far- field data (through Uo n ) are encoded in these scalar 
relationships. However, for problems of size "2+2" and larger, it be- 
comes difficult to make general statements about a particular problem. 
In general, for "2+2" problems, w has elements that are fourth order 
polynomials in \a\. Thus, (|14|) will in general also be a fourth or- 
der polynomials in \a\, with coefficients depending on boundary data 
and global fluxes. We discuss below situations of physical interest 
(including the main example from the last section) where significant 
simplification occurs in these expressions. 
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4.2 Application to physical example 



Let's go back to the concrete example in section [3] where G is given 
in ©. The nullspace of G, Af(G), is: 

/ §r 1 1 \ 
Af(G) T = +1 +1 (16) 

V 00 01110/ 

Each row of the matrix above corresponds to the effect of one global 
flux. We describe the physical interpretations of these global fluxes 
below: 

qi'. This corresponds physically to a heat flux downwards in the lower 
region. This heat flux is balanced by a heat flux in the upper 
region. In addition, the temperature gradient in the lower region 
drives a vapour flux there that is matched by a pressure driven 
gradient flow in the upper region. 

qi'. This corresponds physically to a flux of water downwards in the 
lower region. This flux is generated by condensation at the two 
phase boundary which is provided by vapour flow downwards 
and heat flow upwards in the upper region. 

qy. This term represents an additive shift in temperature and pres- 



The data term b is zero in this example, so the well-posedness is 
determined entirely by the global fluxes. Specifically, the term Uo n in 
(P3|) and (fT5]) is given by 



where by the last term above we mean the span of the column vectors, 
which has rank 2. This is our first simplification: that only two of the 
global fluxes (the first two described above) enter the stability condi- 
tion. Upon reflection, it is clear that shifting the overall temperature 
and pressure (the effect of q^) in this linear problem should not affect 
the stability of the problem. 

We now find the left nullspace w of (GM) T : 
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The first three elements are 0(|a|), corresponding to rows in G rep- 
resenting pure Dirichlet conditions. This is the second simplification 
observed in this example, that the vector w T has relatively simple 
dependence on \a\. 

By simple calculation, we get: 

w G U 0n = 2 \a\ K+ + K- ^ 

to be the well-posedness criteria, that is 

(K+ - K~) qi + (2 + K+ + K~)q 2 + 0. (17) 

The simplifications noted above are what lead to this particularly clear 
form. We note that in the physical situation described in [3], q 2 < 
(liquid moves toward the two phase zone) and q\ « so we expect the 
well-posedness term in (fT7|) to be negative. 

We consider now the choice of residual velocity (fT5|) for this case. 
For instance, if v = e§ (the physical velocity choice), we have: 

. . | (K + — K-) qi + (K+ +K- + 2)q 2 

A = 2|a| K^nr- • 

This is a stable choice for velocity, considering the predicted sign of 
the term (|17() above. However, when we consider a — > oo (finer dis- 
cretizations), we get a very stiff problem and the time step we take 
should be correspondingly small. Therefore, better numerical proper- 
ties are obtained when we choose v such that A is independent of wave 
number. This goal can be achieved by choosing residuals of Dirichlet 
conditions. In this case, we can choose v T = ei (corresponding to the 
saturation condition), then 

{K+ - K-) qi + (K+ +K-+ 2)q 2 
K+ + K- + 2 

and we have an optimal velocity to evolve the interface. Numerical 
validation of the superior performance of this residual velocity over 
the physical one is given in Section [5] below. 

A straight forward calculation shows that for all problems of class 
C identified at the end of section [2] similar results are obtained. That 
is, the vector w in (|14ft and (115p has terms of order a in entries cor- 
responding to pure Dirichlet conditions and constant in entries cor- 
responding to pure Neumann conditions. Also, only the two global 
fluxes corresponding to derivative quantities enter these expressions. 
Similar simplifications are found in other cases when the interface 
conditions separate into pure Dirichlet and pure Neumann conditions. 
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4.3 Extension to the 3-dimensional problem 

The linear analysis is similar except that we assume an initial flat 
interface to be z = 0, and a small perturbation of driving function is 
given: 

GU = b + efe iax e if3y (18) 
then consequently we have updated solution and interface: 

rj ~ 7] + er?i H 

u ~ uo + eui + • • • 

where 

m (x,y)=f ]1 e iax e il3y 

Then the separable solutions of Laplacian equation Au = are in 
form of: 

u = e iax e i(3y u 



where 

u = ce 



and c = (A + , A , B + , B ) T is any arbitrary constant vector deter- 
mined by the fixed boundary conditions. So we have: 

If we expand the solution along the new interface, we have: 

u(x, y, r](x, y)) = u (x, y, 0) + er/iu 02 (x, y, 0) + eui(x, y, 0) + 0(e 2 ) 

(20) 

and since the leading order term of interface is t]q = 0, we can consider 
roughly u n = u z . 

u n (x, y, r](x)) = u 0z (x, y, 0) + e[u lz (x, y, 0) + r/iu 0z2 (x, y, 0)] + 0(e 2 ) 

(21) 

Substitute those expansions into (pOj) . equating coefficients of different 
orders, we have: 

[GUi + me iax e il3y V n} = f e iax e i/3y (22) 

where Uo n = ( u ozz' u 0z) T i s dependent on uo, and Ui = (u[ z ,uJ) T 
with ui in (3.15). 
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The base solution uo in 3-D retains the same form in the 2-D case. 
We need the solutions independent of coordinate variables x and y, 
which lead to Uq = zr 1 + r , and consequently 

(r 1 , r°) T GAA(G) 

we always have Uo n = (0,ri) T . Hence: 

U 0n G S[Af(G)] = ( 3x4 , AA(G) i4 ) T (23) 

where M(G)i4 represents the left four columns of M(G). 

After we figure out Uo n , the next step is to find Ui. We reserve 
the notation of matrix M to provide us a more clear structure of Ui. 

lJ 1 (x,y,0) = (n lz ,u 1 ) = e iax e i ^Mc 

and 

M = / TV^T^h 

where I4 is a 4 x 4 identity matrix. Thus the equation (3.18) becomes 

GMC + TftGUon = f ( 24 ) 

Denoting its left nullspace by w = AA[(GM) r ], we still have the same 
criteria for well-posedness: 

w T GU 0n + (25) 
Applying this to our physical example in 3-D, we have: 

w'G U „ = 2 v^Tffi g - K ' )qi K :f K . + ~ + 2 '" 2 * 

This provides the same information as the 2-D problem. A study of 
the properties of interface velocities also gives results analogous to the 
2-D case. 



5 A 2-D computational example 

In this section, we show how the residual velocities work on the vec- 
tor Laplacian problem by computing the linear example of section [3j 
In addition to the interface conditions, the following fixed boundary 
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conditions are given: On the bottom, we assume the medium is sat- 
urated with liquid water, thus s = 1, and give a constant reference 
temperature 

T-(x,0)=T . 

On the top boundary, we make the boundary impenetrable to vapour 

dP 
on 

and take the heat flux to be given: 

K+^- = h{x). 
on 

We take fi(x) = 2 + sin(x)/2 and To = 10 as our numerical example. 
These conditions lead to a global flux q2 < (water flux upwards 
towards the interface) at each interface location. In the computations 
below, we take K + = K_ so the well-posedness condition (fT7|) is 
satisfied with the sign that makes our velocity sign choices stable. 

5.1 Exact solution 

In order to get an exact solution that can be compared with the nu- 
merical approximation, we start by building two piecewise functions 

f K+T+ in D+ [ P in D + 

U = { V ={ 

\ K-T- - s in D~ [ T~ + s in D~ 

The two Laplacian equations AU = and AV = hold in the whole 
domain discarding the interface. The continuity gives the Neumann 
conditions across the interface: 

K+T+ - K~T~ = -s n and P n = T n + s n 

and because s = at the interface, the Dirichlet conditions across the 
interface are: 

K + T + = K~T~ and P = T~ 

Note that only when K + = the problem is equivalent to our 
original one, so we set K + = K~ and try to find the location of the 
interface s = 0. Solving the problem, we have: 

sin(x) sinh(-y) „ 
T+ = 9 + 2y+ y 7 P = ll 
2 cosh(L) 

sin(a;) sinh(y) sin (a) sinh(y) 

T =W + y + — — — - s = l-y — — - 

4 cosh(L) 4 cosh(L) 
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and the location of s = is given implicitly by 

4(1 — y) cosh(L) — sin(x) sinh(y) = 

By setting x equally distributed in the interval [0, 2ir], we can solve 
this transcendental equation for y, and thus find the coordinates of 
points on the interface to arbitrarily high precision. 

5.2 Finite difference method in mapped rect- 
angular domain 

To handle the two irregular subdomains separated by curved interface, 
we map them into two rectangular domains, so that finite difference 
method can be naturally employed. A finite element method would 
be a good alternative for more irregular domains. A preliminary finite 
element implementation was done in [7] which required considerable 
grid refinement near the interface. The finite difference implementa- 
tion given here demonstrates that this is not necessary for accurate 
solutions. 

Suppose y = h(x) is the known interface, we will do the following 
mapping: 

D + : (x, y) -» (xi = x, Vl = 1 + (y - h(x))/(L - h(x))) 
D~ : (x,y) ->(x 2 =x, y 2 = y/h(x)) 

The new domain is [0, 2tt] x [0, 2], and the interface y = h(x) has been 
mapped into y\ = y 2 = 1. Then by simple calculation, we have: 

3T h'(x 2 )y 2 8T 
dx 2 h(x 2 ) dy 2 

1 dT 
h(x 2 ) dy 2 

Note that = 0, ^7 7^ 0. Rewrite the Laplacian equations in the 



dT _ dT dx 2 dT dy 2 

dx dx 2 dx dy 2 dx 

dT dT dx 2 dT dy 2 

dy dx 2 dy dy 2 dy 
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mapped lower subdomain: 

AT - _ d 2 T d 2 T _ 1 d 2 T dx 2 d ( dT\ h'{x 2 )y 2 dy 2 d ( dT s 



<9y 2 cte 2 h{x2) 2 dy\ dx dx 2 \dx J h(x 2 ) dx dy 2 \dx 
1 d 2 T dx 2 d f dT h'(x 2 )y 2 dT' 



+ 



h(x 2 ) 2 dy\ dx dx 2 \dx 2 h(x 2 ) dy 2 
h'{x 2 )y 2 dy 2 d ( dT h'(x 2 )y 2 dT' 



h(x 2 ) dx dy 2 \dx 2 h(x 2 ) dy 2 
1 + h'(x 2 ) 2 y 2 d 2 T~ d 2 T~ _ y 2 h'(x 2 ) d 2 T~ 



h{x 2 ) 2 dy\ dx\ h(x 2 ) dx 2 dy 2 

| 2h'{x 2 ) 2 - h"{x 2 )h{x 2 ) &T- 
h(x 2 ) 2 m dy 2 

(26) 

Similar calculation gives the Laplacian equation in the mapped upper 
subdomain: 

AT+ = l + h'(x 1 ) 2 (2-y 1 ) 2 d 2 T+ + (PT+_ _ 2 (2- yi )h'( Xl ) d 2 T+ 
(L — h(xi)) 2 dy 2 dx 2 (L — h(x\)) dx\dy\ 

2h'{x l ) 2 + h"{x 1 ){L-h{x 1 )) dT+ 
(L-h( Xl )Y [ Vl) dy x 

(27) 

For the pressure P and saturation s, we get similar formulas. 
The Neumann interface and fixed boundary conditions should be 
reformulated as well in the new coordinates. For the Neumann 
(physical) velocity computation, we choose the residual of mass 
conservation condition P n — T~ — s n = to move the interface. Al- 
ternatively, the Dirichlet condition s = is picked as the Dirichlet 
residual velocity. 

In order to make our numerical approximation second order 
accurate, we use central difference scheme all through the dis- 
cretization. With iV x iV to be the resolution on either side, the 
resultant coefficient matrix is of order 4N 2 x 4N 2 . 



5.3 The numerical results 

We start from iV = 10, and increase it to get higher accuracy. 
L = 2 is fixed all through the computation. To get the results 
in Table [H we set the timestep At = 0.02 for the Neumann 
(physical) velocity, and compute till T = 24, which is long enough 
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Grids 


errlnf 


ratio 


errT 


ratio 


errS 


ratio 


10x10 


0.0031 


*** 


0.0075 


*** 


2.8773e-5 


*** 


20x20 


8.0341e-4 


3.859 


0.0019 


3.947 


6.4194e-6 


4.48 


40x40 


2.0019e-4 


4.013 


4.8504e-4 


3.917 


1.6205e-6 


3.961 



Table 1: Convergence test of Neumann residual velocity 



Grids 


errlnf 


ratio 


errT 


ratio 


errS 


ratio 


10x10 


0.0032 


*** 


0.0075 


*** 


2.9558e-5 


*** 


20x20 


8.0341e-4 


3.983 


0.0019 


3.947 


6.4127e-6 


4.601 


40x40 


2.0019e-4 


4.013 


4.8504e-4 


3.917 


1.6201e-6 


3.960 



Table 2: Convergence test of Dirichlet residual velocity 



to reach the steady state. The errors errlnf, errT, errS are 
respectively the errors of interface, temperature and saturation 
in maximum norm. For the Dirichlet (artificial) residual velocity, 
we take At = 0.2, and iterate till T = 24, the errors are shown 
in Table [21 Note that at (almost) steady state, the errors of the 
two methods are (almost) identical since the discretization of the 
steady-state solution is identical. 

Obviously, both methods have second order accuracy implied 
by the convergent rate. Now we try to demonstrate the advan- 
tage of the Dirichlet (artificial) residual velocity by setting the 
timestep At a free parameter. As we know, since we have to 
use explicit scheme for x n = i?(x) to evolve the interface, the 
timestep has to be reduced when refining the grid to make it 
within the stability region. In Table [31 we observe this prop- 
erty in Neumann residual velocity, but not in Dirichlet one. This 
means, even when we have very fine grids, the timestep can be 
kept reasonably large. Since we know in each time iteration, 
the computational complexity is almost the same, allowing larger 
timesteps can greatly reduce the computational cost. 

In Table [31 VN and VD represent respectively the Neumann 
and Dirichlet residual velocities. It's not hard to see for Neumann 
residual velocity, At = 0.12 is appropriate for 10 x 10 grids to 
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Grids 


timestep 


RV 


errlnf 


errT 


errS 


10x10 


0.20 


VN 


0.9213 


0.8734 


0.1029 


0.15 


VN 


0.0453 


0.0863 


0.0122 


0.12 


VN 


0.0031 


0.0075 


2.8773e-5 


10x10 


0.20 


VD 


0.0032 


0.0075 


2.9558e-5 


20x20 


0.12 


VN 


12.0340 


6.2583 


0.8346 


0.08 


VN 


73.4589 


4.9261 


1.3080 


0.06 


VN 


8.0341e-4 


0.0019 


6.4194e-6 


20x20 


0.20 


VD 


8.0341e-4 


0.0019 


6.4127e-6 


40x40 


0.06 


VN 


9.6403 


24.4088 


0.9594 


0.04 


VN 


3.8939 


3.5538 


0.8762 


0.03 


VN 


2.0019e-4 


4.8504e-4 


1.6205e-6 


40x40 


0.20 


VD 


2.0019e-4 


4.8504e-4 


1.6201e-6 



Table 3: Performance of different residual velocities 



get convergence, and when we refine to 20 x 20, At = 0.06 is 
correspondingly decreased, and for 40 x 40, it becomes At = 0.03. 
However, we can use At = 0.2 for all three resolutions if the 
Dirichlet residual velocity is chosen. Considering the appropriate 
timesteps At = 0.03 and At = 0.2 when the resolution is 40 x 40, 
Dirichlet residual velocity is obviously better than the physical 
Neumann one. 

6 Conclusions and discussion 

Local, linear well-posedness criteria were developed for a general 
class of "2+2" vector Laplacian problems. Much simplification 
can be obtained in certain cases of physical interest. This is 
demonstrated in a simplified model of two-phase flow with phase 
change in porous media. The theory also allows the investigation 
of the stability and numerical properties of residual velocities to 
compute the steady state. It was possible to identify a superior 
velocity choice for the physical example and other problems in 
its general class. The theory was verified computationally using 
a finite difference method with mapped coordinates. In general, 
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this approach offers a way to improve the performance of certain 
steady-state interface computations with very little effort. Ongo- 
ing work suggests that these ideas can be extended to interface 
problems involving more general second order elliptic systems, 
such as problems from viscous, incompressible fluid flow. Ex- 
tension to higher order problems with biharmonic operators are 
also possible. Efficient implementation of the method in a finite 
element framework is still an open problem. 
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